#ifndef SHARED_H_INCLUDED
#define SHARED_H_INCLUDED

using namespace std;

inline double fitness(int fitnessType,double robustness,double complexity)
{
    double fit = 0.0;
    if(fitnessType == 0)
    {
        fit = robustness;
    }
    else if(fitnessType == 1) // select for simplicity
    {
        fit = 1 - complexity; // this is correct. The more Reduced Rules (larger complexity), the more complex this lineage is.
    }
    else if(fitnessType == 2)
    {
        fit = robustness/complexity;
    }
    else if(fitnessType == 3)
    {
        fit = pow(robustness,5)/complexity;
    }
    else if(fitnessType == 4)
    {
        fit = pow(robustness,100)/complexity;
    }
    return 0.000001 + fit;
}

inline string types2Desc(int type, string lineageType)
{
    size_t found;
    char buffer [50];
    string tmpStr;

    found = lineageType.find("exp");
    if (found!=string::npos)
    {
        if(type == Lineages::sparableCellType)
        {
            return "X"; // sparableCellType values will not affect fitness when lost
        }
        else if ( (-1 >= type) and (type >= -8) )
        {
            sprintf(buffer,"%d",(0 - type));
            tmpStr.insert(0,buffer);
            return(tmpStr);
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }

    found = lineageType.find("fun");
    if (found!=string::npos)
    {
        if(type == Lineages::sparableCellType)
        {
            return "X"; // sparableCellType values will not affect fitness when lost
        }
        else if (type == -1)
        {
            return "Neu";
        }
        else if (type == -2)
        {
            return "Bla";// negative values are terminal types
        }
        else if (type == -3)
        {
            return "Epi";
        }
        else if (type == -4)
        {
            return "Ger";
        }
        else if (type == -5)
        {
            return "Gla";
        }
        else if (type == -6)
        {
            return "Int";
        }
        else if (type == -7)
        {
            return "Mus";
        }
        else if (type == -8)
        {
            return "Str";
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }

	found = lineageType.find("cpe"); // "cpe" means "Celegans post embryonic"
    if (found!=string::npos)
    {
        if(type == Lineages::sparableCellType)
        {
            return "X"; // sparableCellType values will not affect fitness when lost
        }
        else if (type == -1)
        {
            return "Neu";
        }
        else if (type == -2)
        {
            return "Bla";// negative values are terminal types
        }
        else if (type == -3)
        {
            return "Epi";
        }
        else if (type == -4)
        {
            return "Ger";
        }
        else if (type == -5)
        {
            return "Gla";
        }
        else if (type == -6)
        {
            return "Int";
        }
        else if (type == -7)
        {
            return "Mus";
        }
        else if (type == -8)
        {
            return "Str";
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }

    found = lineageType.find("cel");
    if (found!=string::npos)
    {
        if(type == Lineages::sparableCellType)
        {
            return "X"; // sparableCellType values will not affect fitness when lost
        }
        else if (type == -1)
        {
            return "ENeu";
        }
        else if (type == -2)
        {
            return "Bla";// negative values are terminal types
        }
        else if (type == -3)
        {
            return "Epi";
        }
        else if (type == -4)
        {
            return "Ger";
        }
        else if (type == -5)
        {
            return "Gla";
        }
        else if (type == -6)
        {
            return "Int";
        }
        else if (type == -7)
        {
            return "Mus";
        }
        else if (type == -8)
        {
            return "Str";
        }
        else if (type == -9)
        {
            return "INeu";
        }
        else if (type == -10)
        {
            return "MNeu";
        }
        else if (type == -11)
        {
            return "PNeu";
        }
        else if (type == -12)
        {
            return "SNeu";
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }

    found = lineageType.find("pma");
    if(found!=string::npos)
    {
        if(type == Lineages::sparableCellType)
        {
            return "X";
        }
        else if(type == -7)
        {
            return "?";
        }
        else if (type == -1)
        {
            return "Ner";
        }
        else if (type ==  -2)
        {
            return "Pha";// negative values are terminal types
        }
        else if (type == -3)
        {
            return "Epi";
        }
        else if (type == -4)
        {
            return "Mus";
        }
        else if (type == -5)
        {
            return "Int";
        }
        else if (type == -6)
        {
            return "Ger";
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }

    found = lineageType.find("hro");
    if(found!=string::npos)
    {
        if(type == -1)
        {
            return "end";
        }
        else if(type == -2)
        {
            return "epi";
        }
        else if (type == -3)
        {
            return "mes";
        }
        else if (type == -4)
        {
            return "mus";
        }
        else if (type == -5)
        {
            return "ner";
        }
        else if (type == -6)
        {
            return "not";
        }
        else if (type == -7)
        {
            return "und";
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }

    found = lineageType.find("mhr");
    if(found!=string::npos)
    {
        if(type == -1)
        {
            return "end";
        }
        else if(type == -2)
        {
            return "epi";
        }
        else if (type == -3)
        {
            return "mes";
        }
        else if (type == -4)
        {
            return "mus";
        }
        else if (type == -5)
        {
            return "ner";
        }
        else if (type == -6)
        {
            return "not";
        }
        else if (type == -7)
        {
            return "und";
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }
    cerr << "An unexpected filename of the lineages!" << endl;
    exit(0);
    return(0);
}

inline int lin2typeCnt(string lineageType)
{
    size_t found;
    found = lineageType.find("fun");
    if(found!=string::npos)
    {
        return 8;
    }
    found = lineageType.find("cpe");
    if(found!=string::npos)
    {
        return 8;
    }
    found = lineageType.find("pma");
    if(found!=string::npos)
    {
        return 7;
    }
    found = lineageType.find("cel");
    if(found!=string::npos)
    {
        return 12;
    }
    found = lineageType.find("hro");
    if(found!=string::npos)
    {
        return 7;
    }
    found = lineageType.find("mhr");
    if(found!=string::npos)
    {
        return 7;
    }
    return(0);
}

inline int desc2Types(string desc, string lineageType)
{
    size_t found;

    found = lineageType.find("exp");
    if (found!=string::npos)
    {
        if(desc == "X")
        {
            return Lineages::sparableCellType; // sparableCellType values will not affect fitness when lost
        }
        else if ( (atoi(desc.c_str() ) >= 1) and (atoi(desc.c_str() ) <= 8) )
        {
            return(0 - atoi(desc.c_str() ));
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }

    found = lineageType.find("fun");
    if (found!=string::npos)
    {
        if(desc == "Dea")
        {
            return Lineages::sparableCellType; // sparableCellType values will not affect fitness when lost
        }
        else if (desc == "X")
        {
            return Lineages::sparableCellType;
        }
        else if (desc ==  "Neu")
        {
            return -1;// negative values are terminal types
        }
        else if (desc == "Bla")
        {
            return -2;
        }
        else if (desc == "Epi")
        {
            return -3;
        }
        else if (desc == "Ger")
        {
            return -4;
        }
        else if (desc == "Gla")
        {
            return -5;
        }
        else if (desc == "Int")
        {
            return -6;
        }
        else if (desc == "Mus")
        {
            return -7;
        }
        else if (desc == "Str")
        {
            return -8;
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }

	found = lineageType.find("cpe");
    if (found!=string::npos)
    {
        if(desc == "Dea")
        {
            return Lineages::sparableCellType; // sparableCellType values will not affect fitness when lost
        }
        else if (desc == "X")
        {
            return Lineages::sparableCellType;
        }
        else if (desc ==  "Neu")
        {
            return -1;// negative values are terminal types
        }
        else if (desc == "Bla")
        {
            return -2;
        }
        else if (desc == "Epi")
        {
            return -3;
        }
        else if (desc == "Ger")
        {
            return -4;
        }
        else if (desc == "Gla")
        {
            return -5;
        }
        else if (desc == "Int")
        {
            return -6;
        }
        else if (desc == "Mus")
        {
            return -7;
        }
        else if (desc == "Str")
        {
            return -8;
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }

    found = lineageType.find("cel");
    if (found!=string::npos)
    {
        if(desc == "Dea")
        {
            return Lineages::sparableCellType; // sparableCellType values will not affect fitness when lost
        }
        else if (desc == "X")
        {
            return Lineages::sparableCellType;
        }
        else if (desc ==  "ENeu")
        {
            return -1;// negative values are terminal types
        }
        else if (desc ==  "INeu")
        {
            return -9;// negative values are terminal types
        }
        else if (desc ==  "MNeu")
        {
            return -10;// negative values are terminal types
        }
        else if (desc ==  "PNeu")
        {
            return -11;// negative values are terminal types
        }
        else if (desc ==  "SNeu")
        {
            return -12;// negative values are terminal types
        }
        else if (desc == "Bla")
        {
            return -2;
        }
        else if (desc == "Epi")
        {
            return -3;
        }
        else if (desc == "Ger")
        {
            return -4;
        }
        else if (desc == "Gla")
        {
            return -5;
        }
        else if (desc == "Int")
        {
            return -6;
        }
        else if (desc == "Mus")
        {
            return -7;
        }
        else if (desc == "Str")
        {
            return -8;
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }

    found = lineageType.find("pma");
    if(found!=string::npos)
    {
        if(desc == "X")
        {
            return Lineages::sparableCellType; // sparableCellType values will not affect fitness when lost
        }
        else if (desc == "?")
        {
            return -7;
        }
        else if (desc ==  "Ner")
        {
            return -1;// negative values are terminal types
        }
        else if (desc == "Pha")
        {
            return -2;
        }
        else if (desc == "Epi")
        {
            return -3;
        }
        else if (desc == "Mus")
        {
            return -4;
        }
        else if (desc == "Int")
        {
            return -5;
        }
        else if (desc == "Ger")
        {
            return -6;
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }

    found = lineageType.find("hro");
    if(found!=string::npos)
    {
        if(desc == "end")
        {
            return -1;
        }
        else if (desc == "epi")
        {
            return -2;
        }
        else if (desc == "mes")
        {
            return -3;
        }
        else if (desc == "mus")
        {
            return -4;
        }
        else if (desc == "ner")
        {
            return -5;
        }
        else if (desc == "not")
        {
            return -6;
        }
        else if (desc == "und")
        {
            return -7;
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }

    found = lineageType.find("mhr");
    if(found!=string::npos)
    {
        if(desc == "end")
        {
            return -1;
        }
        else if (desc == "epi")
        {
            return -2;
        }
        else if (desc == "mes")
        {
            return -3;
        }
        else if (desc == "mus")
        {
            return -4;
        }
        else if (desc == "ner")
        {
            return -5;
        }
        else if (desc == "not")
        {
            return -6;
        }
        else if (desc == "und")
        {
            return -7;
        }
        cerr << "An unexpected cell type is encountered!" << endl;
        exit(0);
    }
    cerr << "An unexpected filename of the lineages!" << endl;
    exit(0);
    return(0);
}

void loadLineageCollection(string file, int skipRows , int loadRows ,vector<vector<int> > & termTypesCollection, vector<vector<string> > & termLineagesCollection)
{
    ifstream ifs;
    string tmpStr;
    ifs.open(file.c_str());
    vector<int> termTypes;
    vector<string> termLineages;
    char c;
    int curRows = 0;
    int loadedRows = 0;
    while(!ifs.eof())
    {
        if(curRows < skipRows)
        {
            ifs.ignore(999999999,'\n');
            curRows ++;
            continue;
        }
        ifs >> tmpStr;
        termTypes.push_back(desc2Types(tmpStr,file));
        ifs.ignore(1000,' ');
        ifs >> tmpStr;
        if(!((tmpStr.substr(0,1) == "0") or (tmpStr.substr(0,1) == "1")))
        {
            cerr << "An unexpected cell lineage description is encountered!" << endl;
            exit(0);
        }
        termLineages.push_back(tmpStr);
        ifs.ignore(1000,';');
        c = ifs.peek();
        if((c == '\n') or (ifs.eof()))
        {
            termTypesCollection.push_back(termTypes);
            termTypes.clear();
            termLineagesCollection.push_back(termLineages);
            termLineages.clear();
            if(c == '\n') ifs.ignore(1,'\n');
            curRows ++;
            loadedRows ++;
            if(loadedRows >= loadRows) break;
        }
    }
    ifs.close();

    return;
}

void writeLineageCollection(string file,vector<vector<int> > & termTypesCollection, vector<vector<string> > & termLineagesCollection)
{
    ofstream ofs;
    ofs.open(file.c_str(),ios_base::out);
    vector<int> termTypes;
    vector<string> termLineages;
    for(unsigned int i=0;i<termTypesCollection.size();i++)
    {
        for(unsigned int j=0;j<termTypesCollection[i].size();j++)
        {
            ofs << types2Desc(termTypesCollection[i][j],file);
            ofs << ' ';
            ofs << termLineagesCollection[i][j];
            ofs << ' ';
            ofs << ';';
        }
        ofs << '\n';
    }

    ofs.close();

    return;
}

void appendLineageCollection(string file,vector<vector<int> > & termTypesCollection, vector<vector<string> > & termLineagesCollection)
{
    ofstream ofs;
    ofs.open(file.c_str(),ios_base::app);
    vector<int> termTypes;
    vector<string> termLineages;
    for(unsigned int i=0;i<termTypesCollection.size();i++)
    {
        for(unsigned int j=0;j<termTypesCollection[i].size();j++)
        {
            ofs << types2Desc(termTypesCollection[i][j],file);
            ofs << ' ';
            ofs << termLineagesCollection[i][j];
            ofs << ' ';
            ofs << ';';
        }
        ofs << '\n';
    }

    ofs.close();

    return;
}

inline double random(double start, double end)
{
    return start +(end-start)*rand()/(RAND_MAX + 1.0);
}

void loadRealLineages(string file,vector<int> & terTypes,vector<string> & terLineages)
{
    ifstream ifs;
    string tmpStr;
    ifs.open(file.c_str());
    ifs.ignore(1000,'\n'); // ignore the header
    /*
    Lineage	Name	Class
    0000000000	AINL	Neu
    0000000001	Dead	Dea
    */
    while(ifs >> tmpStr)
    {
    	//cerr << tmpStr << endl;
		if(tmpStr.substr(0,1) == "#") {
			ifs.ignore(1000,'\n');
			continue;
		}
        if(!((tmpStr.substr(0,1) == "0") or (tmpStr.substr(0,1) == "1")))
        {
            cerr << "An unexpected cell lineage description is encountered!" << endl;
            exit(0);
        }
        terLineages.push_back(tmpStr);
        ifs >> tmpStr;
        ifs >> tmpStr;
        terTypes.push_back(desc2Types(tmpStr,file));
        ifs.ignore(1000,'\n');
    }
    ifs.close();

    return;
}

bool startWith(string str1, string str2) // is str1 start with str2?
{
    string tmpStr = str1.substr(0,(int)str2.length());
    if(tmpStr == str2) return true;
    else return false;
}

void outputLineages(string file,const vector<int> cellTypes,const vector<string> cellLineages)
{
    ofstream ofs;
    ofs.open(file.c_str());
    ofs << "Lineage\tName\tClass" << endl;
    for(unsigned int i=0; i<cellLineages.size(); i++)
    {
        ofs << cellLineages[i] << "\t" << cellTypes[i] << "\t" << cellTypes[i] << "\n";
    }
    ofs.close();
    return;
}

void generateRandomLineageTopologyLimitDepth(int numTerm, int depthLimit, vector<string> & res)
{
    if(pow(2,depthLimit) < numTerm)
    {
        cerr << "Depth "<< depthLimit << " could not contain ";
        cerr << numTerm << " terminal cells" << endl;
        exit(3);
    }
    int tmpTerm,maxDepth;
    string newTerm1,newTerm2;

    maxDepth = 0;

    while(maxDepth != depthLimit)
    {
        res.clear();
        res.push_back("");
        for(int i=1; i<numTerm; i++)
        {
            tmpTerm = (int)random(0,res.size());
            newTerm1 = res[tmpTerm];
            if(newTerm1.length() >= (unsigned)depthLimit)
            {
                i --;
                continue;
            }
            newTerm1 += "1";
            res.push_back(newTerm1);
            newTerm2 = res[tmpTerm];
            newTerm2 += "0";
            res[tmpTerm] = newTerm2;
            if(maxDepth < (signed)newTerm1.length())
            {
                maxDepth = newTerm1.length();
            }
        }
    }
    random_shuffle(res.begin(),res.end());
    return;
}

void generateRandomLineageTopology(int numTerm, vector<string> & res)
{
    int tmpTerm;
    string newTerm1,newTerm2;

    res.clear();
    res.push_back("");
    for(int i=1; i<numTerm; i++)
    {
        tmpTerm = (int)random(0,res.size());
        newTerm1 = res[tmpTerm];
        newTerm1 += "1";
        res.push_back(newTerm1);
        newTerm2 = res[tmpTerm];
        newTerm2 += "0";
        res[tmpTerm] = newTerm2;
    }
    random_shuffle(res.begin(),res.end());
    return;
}

void depthFixedShuffle(const vector<string> & termLineages, vector<string> & resLineages)
{
    unsigned int maxDepth,minDepth;
    map<int, vector<int> > depth2Term;
    vector<int> tmpVec;
    maxDepth = 0;
    minDepth = 999;
    resLineages.clear();

    for(unsigned int i=0;i<termLineages.size();i++) {
        if(termLineages[i].length() > maxDepth) maxDepth = termLineages[i].length();
        if(termLineages[i].length() < minDepth) minDepth = termLineages[i].length();
        depth2Term[termLineages[i].length()].push_back(i);
        resLineages.push_back(termLineages[i]);
    }

    for(unsigned int i=minDepth;i<=maxDepth;i++) {
        tmpVec = depth2Term[i];
        random_shuffle(tmpVec.begin(),tmpVec.end());
        for(unsigned int j=0;j<depth2Term[i].size();j++)
        {
            resLineages[depth2Term[i][j]] = termLineages[tmpVec[j]];
        }
    }

    return;
}

void depthSemifixedShuffle(const vector<string> & termLineages, vector<string> & resLineages, double probOutOfDepth)
{
    unsigned int maxDepth,minDepth,tmpDepth1,tmpDepth2,tmpInt1,tmpInt2,tmpSwap;
    map<int, vector<int> > depth2Term;
    vector<int> tmpVec;
    unsigned int numOutOfDepth = int(termLineages.size() * probOutOfDepth);
    maxDepth = 0;
    minDepth = 999;
    resLineages.clear();

    for(unsigned int i=0;i<termLineages.size();i++) {
        if(termLineages[i].length() > maxDepth) maxDepth = termLineages[i].length();
        if(termLineages[i].length() < minDepth) minDepth = termLineages[i].length();
        depth2Term[termLineages[i].length()].push_back(i);
        resLineages.push_back(termLineages[i]);
    }

	for(unsigned int i=0;i<numOutOfDepth;i++) {
		tmpDepth1 = (tmpDepth2 = 5);
		while(tmpDepth1 == tmpDepth2) {
			tmpDepth1 = (int)random(minDepth,maxDepth);
			tmpDepth2 = (int)random(minDepth,maxDepth);
		}
		tmpVec = depth2Term[tmpDepth1];
		if(tmpVec.size() == 0) {continue;}
		tmpInt1 = int(random((double)0,double(tmpVec.size())));
		tmpVec = depth2Term[tmpDepth2];
		if(tmpVec.size() == 0) {continue;}
		tmpInt2 = int(random((double)0,double(tmpVec.size())));
		tmpSwap = depth2Term[tmpDepth1][tmpInt1];
		depth2Term[tmpDepth1][tmpInt1] = depth2Term[tmpDepth2][tmpInt2];
		depth2Term[tmpDepth2][tmpInt2] = tmpSwap;
	}

    for(unsigned int i=minDepth;i<=maxDepth;i++) {
        tmpVec = depth2Term[i];
        random_shuffle(tmpVec.begin(),tmpVec.end());
        for(unsigned int j=0;j<depth2Term[i].size();j++)
        {
            resLineages[depth2Term[i][j]] = termLineages[tmpVec[j]];
        }
    }

    return;
}

void depthFixedShuffle_twin(const vector<string> & termLineages, const vector<int> termTypes, vector<string> & resLineages)
{
	// within each depth
	// build a twin(twin means same type sibling here) list
	// shuffle the term-types
    unsigned int maxDepth,minDepth;
    map<int, vector<int> > depth2term;
    map<int, int> term2twin; // mapping the terminal id to its twin id
    map<int, int>::iterator it2,it3;
    vector<int> tmpVec,tmpVec2;
    maxDepth = 0;
    minDepth = 999;
    resLineages.clear();
    string tmpTwin1,tmpTwin2,tmpParent;
    map<string,int> lineage2id;
    map<string,int>::iterator it;
    vector<int> unfilledSlot;
    vector<int>::iterator it4;
    int tmpInt;

    for(unsigned int i=0;i<termLineages.size();i++) {
        if(termLineages[i].length() > maxDepth) maxDepth = termLineages[i].length();
        if(termLineages[i].length() < minDepth) minDepth = termLineages[i].length();
        depth2term[termLineages[i].length()].push_back(i);
        resLineages.push_back(termLineages[i]);
        lineage2id[termLineages[i]] = i;

        // populate term2Twin
		tmpTwin1 = termLineages[i];
		tmpParent = tmpTwin1.substr(0,int(tmpTwin1.length()-1));
		tmpTwin2 = tmpTwin1.substr(int(tmpParent.length()),1);
		if(tmpTwin2 == "1")
		{
			tmpTwin2 = "0";
			tmpTwin2.insert(0,tmpParent);
		}
		else
		{
			tmpTwin2 = "1";
			tmpTwin2.insert(0,tmpParent);
		}
		it = lineage2id.find(tmpTwin2);
		if(it == lineage2id.end()) continue; // could not find its twin
		term2twin[i] = lineage2id[tmpTwin2];
		term2twin[lineage2id[tmpTwin2]] = i;
    }

    for(unsigned int i=minDepth;i<=maxDepth;i++) {
        tmpVec = depth2term[i];
        unfilledSlot = tmpVec;
        random_shuffle(unfilledSlot.begin(),unfilledSlot.end());
        random_shuffle(tmpVec.begin(),tmpVec.end());
        // put those with twin at the beginning.
        tmpVec2.clear();
		for(it4=tmpVec.begin();it4!=tmpVec.end();it4++) {
			tmpInt = *it4;
			it2 = term2twin.find(tmpInt);
			if((it2 != term2twin.end()) && (termTypes[it2->first] == termTypes[it2->second]))
			{
				tmpVec2.insert(tmpVec2.begin(),tmpInt);
			}
			else
			{
				tmpVec2.push_back(tmpInt);
			}
		}
		tmpVec = tmpVec2;
        while(tmpVec.size() > 0) { // iterate the lineal terminals
			it2 = term2twin.find(tmpVec[0]); // does this terminals have twin terminal?
			if((it2 == term2twin.end()) || (termTypes[it2->first] != termTypes[it2->second]))
        	{   // no, just put it to another no-twin terminals
        		// can do this only because I've handle the twins already
				resLineages[tmpVec[0]] = termLineages[unfilledSlot[0]];
				tmpVec.erase(tmpVec.begin());
				unfilledSlot.erase(unfilledSlot.begin());
        	}
        	else // yes, put it and its twin into another twin position
        	{
				for(it4=unfilledSlot.begin();it4!=unfilledSlot.end();it4++)
				{
					it3 = term2twin.find(*it4);
					if(it3 != term2twin.end())
					{
						resLineages[tmpVec[0]] = termLineages[*it4];
						resLineages[it2->second] = termLineages[it3->second];
						unfilledSlot.erase(it4);
						for(it4=unfilledSlot.begin();it4!=unfilledSlot.end();it4++)
						{
							if(*it4 == it3->second)
							{
								unfilledSlot.erase(it4);
								break;
							}
						}
						tmpVec.erase(tmpVec.begin());
						for(it4=tmpVec.begin();it4!=tmpVec.end();it4++)
						{
							if(*it4 == it2->second)
							{
								tmpVec.erase(it4);
								break;
							}
						}
						break;
					}
				}
        	}
        }
    }

    return;
}

void depthFixedCoalescent(const vector<string> & termLineages,vector<string> & resLineages)
{
    int maxDepth,currentDepth;
    vector<int> nodesInCurrentDepth;
    vector<int> parentsOfPreviousDepth;
    vector<vector<int> > descendantsOfCurrentNode;
    vector<vector<int> > descendantsOfPreviousNode;
    vector<int> tmpVecInt;
    vector<int> aggreDescendant;
    vector<int> withinDepthIndex;
    int tmpCnt = termLineages.size();
    string tmpStr;

    resLineages.clear();

    maxDepth = 0;
    resLineages.clear();
    for(unsigned int i = 0; i < termLineages.size(); i++)
    {
        if((signed)termLineages[i].length() > maxDepth) maxDepth = (signed)termLineages[i].length();
        resLineages.push_back("");
    }

    for(currentDepth = maxDepth; currentDepth > 0; currentDepth --)
    {
        nodesInCurrentDepth.clear();
        nodesInCurrentDepth.insert(nodesInCurrentDepth.begin(),parentsOfPreviousDepth.begin(),parentsOfPreviousDepth.end());
        for(unsigned int i = 0; i<termLineages.size(); i++)
        {
            if(termLineages[i].length() == (unsigned)currentDepth) nodesInCurrentDepth.push_back(i);
        }

        descendantsOfCurrentNode.clear();
        descendantsOfCurrentNode.insert(descendantsOfCurrentNode.begin(),descendantsOfPreviousNode.begin(),descendantsOfPreviousNode.end());
        descendantsOfPreviousNode.clear();
        parentsOfPreviousDepth.clear();
        withinDepthIndex.clear();
        for(unsigned int i=0;i<nodesInCurrentDepth.size();i++)
        {
            withinDepthIndex.push_back(i);
        }

        // randomly picked two
        // record their descendants
        // change lineage string of their descendants
        // make their parent
        // aggregate their descendant for their parent
        random_shuffle(withinDepthIndex.begin(),withinDepthIndex.end());
        for(unsigned int i=0; i<(withinDepthIndex.size()/2); i++)
        {
            aggreDescendant.clear();
            for(unsigned int j=i*2; j<=i*2+1; j++)
            {
                tmpVecInt.clear();
                if(nodesInCurrentDepth[withinDepthIndex[j]] < (signed)termLineages.size())
                {
                    // if they are terminal
                    tmpVecInt.push_back(nodesInCurrentDepth[withinDepthIndex[j]]);
                }
                else
                {
                    //if they are internal
                    tmpVecInt.insert(tmpVecInt.begin(),descendantsOfCurrentNode[withinDepthIndex[j]].begin(),descendantsOfCurrentNode[withinDepthIndex[j]].end());
                }
                for(unsigned int k = 0; k<tmpVecInt.size(); k++)
                {
                    std::stringstream ss;
                    ss << (j - i*2);
                    resLineages[tmpVecInt[k]].insert(0,ss.str());
                    ss.clear();
                }
                aggreDescendant.insert(aggreDescendant.end(),tmpVecInt.begin(),tmpVecInt.end());
            }
            parentsOfPreviousDepth.push_back(tmpCnt);
            descendantsOfPreviousNode.push_back(aggreDescendant);

            tmpCnt ++;
        }
    }
    return;
}

void depthFixedCoalescent2(const vector<int> & termDepths,vector<string> & resLineages)
{
    int maxDepth,currentDepth;
    vector<int> nodesInCurrentDepth;
    vector<int> parentsOfPreviousDepth;
    vector<vector<int> > descendantsOfCurrentNode;
    vector<vector<int> > descendantsOfPreviousNode;
    vector<int> tmpVecInt;
    vector<int> aggreDescendant;
    vector<int> withinDepthIndex;
    int tmpCnt = termDepths.size();
    string tmpStr;

    resLineages.clear();

    maxDepth = 0;
    resLineages.clear();
    for(unsigned int i = 0; i < termDepths.size(); i++)
    {
        if(termDepths[i] > maxDepth) maxDepth = termDepths[i];
        resLineages.push_back("");
    }

    for(currentDepth = maxDepth; currentDepth > 0; currentDepth --)
    {
        nodesInCurrentDepth.clear();
        nodesInCurrentDepth.insert(nodesInCurrentDepth.begin(),parentsOfPreviousDepth.begin(),parentsOfPreviousDepth.end());
        for(unsigned int i = 0; i<termDepths.size(); i++)
        {
            if(termDepths[i] == currentDepth) nodesInCurrentDepth.push_back(i);
        }

        descendantsOfCurrentNode.clear();
        descendantsOfCurrentNode.insert(descendantsOfCurrentNode.begin(),descendantsOfPreviousNode.begin(),descendantsOfPreviousNode.end());
        descendantsOfPreviousNode.clear();
        parentsOfPreviousDepth.clear();
        withinDepthIndex.clear();
        for(unsigned int i=0;i<nodesInCurrentDepth.size();i++)
        {
            withinDepthIndex.push_back(i);
        }

        // randomly picked two
        // record their descendants
        // change lineage string of their descendants
        // make their parent
        // aggregate their descendant for their parent
        random_shuffle(withinDepthIndex.begin(),withinDepthIndex.end());
        for(unsigned int i=0; i<(withinDepthIndex.size()/2); i++)
        {
            aggreDescendant.clear();
            for(unsigned int j=i*2; j<=i*2+1; j++)
            {
                tmpVecInt.clear();
                if(nodesInCurrentDepth[withinDepthIndex[j]] < (signed)termDepths.size())
                {
                    // if they are terminal
                    tmpVecInt.push_back(nodesInCurrentDepth[withinDepthIndex[j]]);
                }
                else
                {
                    //if they are internal
                    tmpVecInt.insert(tmpVecInt.begin(),descendantsOfCurrentNode[withinDepthIndex[j]].begin(),descendantsOfCurrentNode[withinDepthIndex[j]].end());
                }
                for(unsigned int k = 0; k<tmpVecInt.size(); k++)
                {
                    std::stringstream ss;
                    ss << (j - i*2);
                    resLineages[tmpVecInt[k]].insert(0,ss.str());
                    ss.clear();
                }
                aggreDescendant.insert(aggreDescendant.end(),tmpVecInt.begin(),tmpVecInt.end());
            }
            parentsOfPreviousDepth.push_back(tmpCnt);
            descendantsOfPreviousNode.push_back(aggreDescendant);

            tmpCnt ++;
        }
    }
    return;
}

void demirror(vector<string> &fromLin, vector<int> &fromTerm, vector<string> &resLin, vector<int> &resTerm)
// this function just extract a subtree, whose lineage string starts with "0"
// It does not check whether the lineage is symmetric or not
{
    string subLinPrefix = "0";
    resTerm.clear();
    resLin.clear();

    for(unsigned int i=0;i<fromLin.size();i++)
    {
        if(fromLin[i].substr(0,1) == subLinPrefix)
        {
            resLin.push_back(fromLin[i].substr(1));
            resTerm.push_back(fromTerm[i]);
        }
    }
    return;
}

void mirror(vector<string> &fromLin, vector<int> &fromTerm, vector<string> &resLin, vector<int> &resTerm)
// This function just mirror the lineage to get a new lineage with doubled size
{
    string tmpStr;

    resLin.clear();
    resTerm.clear();

    for(unsigned int i=0;i<fromLin.size();i++)
    {
        tmpStr = "0";
        tmpStr += fromLin[i];
        resLin.push_back(tmpStr);
        resTerm.push_back(fromTerm[i]);

        tmpStr = "1";
        tmpStr += fromLin[i];
        resLin.push_back(tmpStr);
        resTerm.push_back(fromTerm[i]);
    }

    return;
}

bool goodDepthSet(const vector<int> & depths) // could this build a tree?
{
	map<int,int> depth2Count;
	map<int,int>::iterator it;
	int tmpInt,minDepth,maxDepth,internCount;
	maxDepth = 0;
	minDepth = 999;
	for(unsigned int i=0;i<depths.size();i++)
	{
		tmpInt = depths[i];
		if(tmpInt > maxDepth) maxDepth = tmpInt;
		if(tmpInt < minDepth) minDepth = tmpInt;
		it = depth2Count.find(tmpInt);
		if(it != depth2Count.end() ) // found
		{
			depth2Count[tmpInt] ++;
		}
		else
		{
			depth2Count[tmpInt] = 1;
		}
	}

	// could this set of depth form a tree?
	internCount = 0;
	for(int i = maxDepth;i>0;i--)
	{
		it = depth2Count.find(i);
		tmpInt = it->second;
		tmpInt += internCount;
		if(tmpInt % 2 == 0 ) // ok
		{
			internCount = tmpInt / 2;
		}
		else
		{
			return false;
		}
	}
	return true;
}

void depthFixedClusterShuffle(const vector<string> & termLineages, const vector<int> & termTypes, const int & biasType, const double & biasPercent, vector<string> & resLineages, vector<int> & resTypes)
// a fraction ("biasPercent") of a certain type("biasType") of cells would be put as clustered as possible
{
    unsigned int maxDepth,minDepth;
    map<int, vector<int> > depth2Term;
    map<int, int> nInTypePerDepth;
    vector<int> tmpVec;
    vector<int> tmpResTypes;
    int nArraged, needArrage, tmpInt;

    maxDepth = 0;
    minDepth = 999;
    for(unsigned int i=0;i<termLineages.size();i++) resLineages.push_back(termLineages[i]);
    sort(resLineages.begin(),resLineages.end());
    tmpResTypes.resize(termTypes.size(),0);
    resTypes.resize(termTypes.size(),0);
    for(unsigned int i=0;i<resLineages.size();i++) {
        for(unsigned int j=0;j<termLineages.size();j++) {
            if(resLineages[i] == termLineages[j]) {
                tmpResTypes[i] = termTypes[j];
            }
        }
    }

    for(unsigned int i=0;i<resLineages.size();i++) {
        tmpInt = resLineages[i].length(); // depth of this cell
        if(tmpInt > (signed)maxDepth) {
            maxDepth = tmpInt;
            nInTypePerDepth[tmpInt] = 0;
        }
        if(tmpInt < (signed)minDepth) {
            minDepth = tmpInt;
            nInTypePerDepth[tmpInt] = 0;
        }
        depth2Term[tmpInt].push_back(i);
        if(tmpResTypes[i] == biasType) nInTypePerDepth[tmpInt] ++;
    }

    for(unsigned int i=minDepth;i<=maxDepth;i++) {
        tmpVec = depth2Term[i];
        random_shuffle(tmpVec.begin(),tmpVec.end());
        needArrage = int(nInTypePerDepth[i] * biasPercent);
        nArraged = 0;
        tmpInt = 0;
        while(nArraged < needArrage) {
            if(tmpResTypes[tmpVec[tmpInt]] == biasType) {
                tmpVec.insert(tmpVec.begin(),tmpVec[tmpInt]);
                tmpVec.erase(tmpVec.begin() + tmpInt + 1);
                nArraged ++;
            }
            tmpInt ++;
        }
        for(unsigned int j=0;j<depth2Term[i].size();j++)
        {
            resTypes[depth2Term[i][j]] = tmpResTypes[tmpVec[j]];
        }
    }

    return;
}

void depthFixedAllClusterShuffle(const vector<string> & termLineages, const vector<int> & termTypes, const double & biasPercent, vector<string> & resLineages, vector<int> & resTypes)
// a fraction ("biasPercent") of all types of cells would be put as clustered as possible
{
    unsigned int maxDepth,minDepth;
    map<int, vector<int> > depth2Term;
    map<int, vector<int> > nInTypePerDepth;
    map<int, int> allTypes;
    vector<int> sortedTypes;
    vector<int> tmpVec;
    vector<int> tmpResTypes;
    int needArrange, tmpInt;

    maxDepth = 0;
    minDepth = 999;
    for(unsigned int i=0;i<termLineages.size();i++) resLineages.push_back(termLineages[i]);
    sort(resLineages.begin(),resLineages.end());
    tmpResTypes.resize(termTypes.size(),0);
    resTypes.resize(termTypes.size(),0);
    for(unsigned int i=0;i<resLineages.size();i++) {
    	allTypes[termTypes[i]] = 1;
        for(unsigned int j=0;j<termLineages.size();j++) {
            if(resLineages[i] == termLineages[j]) {
                tmpResTypes[i] = termTypes[j];
                // Make sure the lineages are sorting alphabetically
                // and the types are sorted accordingly
            }
        }
    }
    for(map<int, int>::iterator it=allTypes.begin();it!=allTypes.end();it++ )
    {
		sortedTypes.push_back(it->first);
    }
    random_shuffle(sortedTypes.begin(),sortedTypes.end());
    // types should be sorted according to sortedTypes in every depth

    for(unsigned int i=0;i<resLineages.size();i++) {
        tmpInt = resLineages[i].length(); // depth of this cell
        if(tmpInt > (signed)maxDepth) {
            maxDepth = tmpInt;
            for(unsigned int j=0;j<sortedTypes.size();j++)
            {
				nInTypePerDepth[sortedTypes[j]].resize(maxDepth+1,0);
            }
        }
        if(tmpInt < (signed)minDepth) {
            minDepth = tmpInt;
        }
        depth2Term[tmpInt].push_back(i);
        tmpVec.clear();
        tmpVec = nInTypePerDepth[tmpResTypes[i]];
        nInTypePerDepth[tmpResTypes[i]][tmpInt] += 1;
    } // how many cells in each type in each depth

    for(unsigned int i=minDepth;i<=maxDepth;i++)
    {
    	tmpVec.clear();
		for(unsigned int k=0;k<sortedTypes.size();k++)
		{
			if(sortedTypes[k] == Lineages::sparableCellType) continue;
			needArrange = int(nInTypePerDepth[sortedTypes[k]][i] * biasPercent);
			nInTypePerDepth[sortedTypes[k]][i] -= needArrange;
			for(int l=0;l<needArrange;l++)
			{
				tmpVec.push_back(sortedTypes[k]);
			}
		}
		for(unsigned int k=0;k<sortedTypes.size();k++)
		{
			needArrange = nInTypePerDepth[sortedTypes[k]][i];
			for(int l=0;l<needArrange;l++)
			{
				tmpVec.insert(tmpVec.begin() + random(0,tmpVec.size()), sortedTypes[k]);
			}
		}
		for(unsigned int j=0;j<depth2Term[i].size();j++)
		{
			resTypes[depth2Term[i][j]] = tmpVec[j];
		}
    }

    return;
}

#endif // SHARED_H_INCLUDED
